home
***
CD-ROM
|
disk
|
FTP
|
other
***
search
/
PD Collection CD 1
/
PD Collection CD 1.iso
/
programer2
/
pari2
/
pari
/
c
/
bibli2
< prev
next >
Wrap
Text File
|
1991-11-28
|
23KB
|
895 lines
/********************************************************************/
/********************************************************************/
/** **/
/** BIBLIOTHEQUE MATHEMATIQUE **/
/** (deuxieme partie) **/
/** **/
/** Copyright Babe Cool **/
/** **/
/********************************************************************/
/********************************************************************/
#include "genpari.h"
/********************************************************************/
/********************************************************************/
/** **/
/** ITERATIONS **/
/** **/
/********************************************************************/
/********************************************************************/
GEN forpari(ep,a,b,ch)
entree *ep;
GEN a,b;
char *ch;
{
GEN p1;
long av=avma,tx=typ(a),la,lx,i,e;
if(tx>2) err(forer1);
if(tx==1)
{
la=lx=lgef(a);if(lx<3) lx=3;
if(!gcmp0(b)) {e=(gexpo(b)>>5)+3;if(e>lx) lx=e;}
p1=newbloc(lx);for(i=0;i<la;i++) p1[i]=a[i];
setlg(p1,lx);p1[-1]=(long)ep->value;ep->value=(void*)p1;
}
else newvalue(ep,a);
p1=(GEN)ep->value;
while (gcmp(p1,b)<=0)
{
lisseq(ch);gaddgsz(p1,1,p1);avma=av;
}
killvalue(ep);return gnil;
}
GEN forstep(ep,a,b,s,ch)
entree *ep;
GEN a,b,s;
char *ch;
{
GEN p1;
long av=avma,tx=typ(a),lx,la,s1,i,e;
s1=gsigne(s);if(!s1) err(forer2);
if(tx>2) err(forer1);
if(tx==1)
{
la=lx=lgef(a);if(lx<3) lx=3;
if(!gcmp0(b)) {e=(gexpo(b)>>5)+3;if(e>lx) lx=e;}
if(typ(s)==2) {affir(a,p1=cgetr(lg(s)));newvalue(ep,p1);avma=av;}
else
{
p1=newbloc(lx);for(i=0;i<la;i++) p1[i]=a[i];
setlg(p1,lx);p1[-1]=(long)ep->value;ep->value=(void*)p1;
}
}
else newvalue(ep,a);
p1=(GEN)ep->value;
if(s1>0)
while (gcmp(p1,b)<= 0)
{
lisseq(ch);gaddz(p1,s,p1);avma = av;
}
else
while (gcmp(p1,b)>= 0)
{
lisseq(ch);gaddz(p1,s,p1);avma = av;
}
killvalue(ep);return gnil;
}
GEN forprime(ep,a,b,ch)
entree *ep;
GEN a,b;
char *ch;
{
GEN p1;
long av=avma,prime=0;
byteptr p=diffptr;
newvalue(ep,gun); p1 = (GEN)ep->value;
while((*p)&&gcmpgs(a,prime)>0) prime += *p++;
if(!*p) err(recprimer);
affsi(prime,p1);
while(gcmp(p1,b)<=0)
{
if(!*p) err(recprimer);
lisseq(ch);addsiz(*p++,p1,p1);avma=av;
}
killvalue(ep);return gnil;
}
GEN fordiv(ep,a,ch)
entree *ep;
GEN a;
char *ch;
{
long i,l,av2,av=avma;
GEN p1,t=divisors(a);
l=lg(t);
newvalue(ep,a); p1=(GEN)ep->value;
av2=avma;
for(i=1;i<l;i++)
{
gaffect(t[i],p1);
lisseq(ch);
avma=av2;
}
avma=av;
killvalue(ep);return gnil;
}
/********************************************************************/
/********************************************************************/
/** **/
/** SOMMES,PRODUITS,VECTEURS,MATRICES ET RECURRENCES **/
/** **/
/********************************************************************/
/********************************************************************/
GEN somme(ep,x,a,b,ch)
entree *ep;
GEN x,a,b;
char *ch;
{
GEN y,z,p1;
long av = avma,tetpil,limite=(av+bot)/2;
newvalue(ep,gun); p1 = (GEN)ep->value; gaffect(a, p1);
y = x;
tetpil = avma;
if(gcmp(a,b)>0)
{
if(!gcmp1(gsub(a,b))) err(sommeer1);
avma = av;
y=gcopy(x);
}
else
do
{
if(avma<limite) {tetpil=avma; y=gerepile(av,tetpil,gcopy(y));}
z=lisexpr(ch);
tetpil=avma; y=gadd(y,z);
gaddgsz(p1,1,p1);
}
while(gcmp(p1,b)<=0);
killvalue(ep);
return gerepile(av,tetpil,y);
}
GEN suminf(ep,a,ch,prec)
entree *ep;
GEN a;
char *ch;
long prec;
{
GEN y,z,p1;
long av=avma,tetpil,limite=(bot+avma)/2,fl=0;
newvalue(ep,gun); p1 = (GEN)ep->value; gaffect(a, p1);
affsr(1,y=cgetr(prec));
do
{
if (avma < limite) {tetpil = avma; y=gerepile(av,tetpil,gcopy(y));}
z=lisexpr(ch); y=gadd(y,z);
gaddgsz(p1,1,p1);
if((!gcmp0(z))&&(gexpo(z)>gexpo(y)-32*(prec-2)-5)) fl=0;
else fl++;
}
while(fl<3);
killvalue(ep);
tetpil=avma; return gerepile(av,tetpil,gsubgs(y,1));
}
GEN produit(ep,x,a,b,ch)
entree *ep;
GEN x,a,b;
char *ch;
{
GEN y,z,p1;
long av=avma,tetpil,limite=(av+bot)/2;
newvalue(ep,gun); p1 = (GEN)ep->value; gaffect(a,p1);
y = x;
tetpil = avma;
if(gcmp(a,b)>0)
{
if(!gcmp1(gsub(a,b))) err(proder1);
avma=av;y=gcopy(x);
}
else
do
{
if(avma<limite) {tetpil=avma; y=gerepile(av,tetpil,gcopy(y));}
z=lisexpr(ch);
tetpil=avma;y=gmul(y,z);
gaddgsz(p1,1,p1);
}
while(gcmp(p1,b)<=0);
killvalue(ep);
return gerepile(av,tetpil,y);
}
GEN prodinf(ep,a,ch,prec)
entree *ep;
GEN a;
char *ch;
long prec;
{
GEN y,z,p1;
long av=avma,tetpil,limite=(av+bot)/2,fl=0;
newvalue(ep, gun); p1=(GEN)ep->value; gaffect(a,p1);
affsr(1,y=cgetr(prec));
do
{
if(avma<limite) {tetpil=avma; y=gerepile(av,tetpil,gcopy(y));}
z=lisexpr(ch); y=gmul(y,z);
gaddgsz(p1,1,p1);
if(gexpo(gsubgs(z,1))>-32*(prec-2)-5) fl=0;else fl++;
}
while(fl<3);
killvalue(ep);
tetpil = avma; return gerepile(av,tetpil,gcopy(y));
}
GEN prodinf1(ep,a,ch,prec)
entree *ep;
GEN a;
char *ch;
long prec;
{
GEN y,z,p1,p2;
long av=avma,tetpil,limite=(av+bot)/2,fl=0;
newvalue(ep,gun); p1=(GEN)ep->value; gaffect(a,p1);
affsr(1,y=cgetr(prec));
do
{
if(avma<limite) {tetpil=avma; y=gerepile(av,tetpil,gcopy(y));}
p2=lisexpr(ch);z=gaddsg(1,p2);
y=gmul(y,z);
gaddgsz(p1,1,p1);
if((!gcmp0(z))&&(gexpo(p2)>-32*(prec-2)-5)) fl=0;else fl++;
}
while(fl<3);
killvalue(ep);
tetpil = avma; return gerepile(av,tetpil,gcopy(y));
}
GEN prodeuler(ep,a,b,ch,prec)
entree *ep;
GEN a,b;
char *ch;
long prec;
{
GEN y,z,p1;
long av=avma,tetpil,prime=0;
byteptr p=diffptr;
newvalue(ep,gun); p1 = (GEN)ep->value;
while((*p)&&gcmpgs(a,prime)>0) prime += *p++;
if(!*p) err(recprimer);
affsi(prime,p1);affsr(1,y=cgetr(prec));
if(gcmp(p1,b)>0)
{
if(!gcmp1(gsub(a,b))) err(recer1);
tetpil=avma;
y=gcopy(y);
}
else
do
{
if(!*p) err(recprimer);
z=lisexpr(ch);
tetpil=avma; y=gmul(y,z);
addsiz(*p++,p1,p1);
}
while(gcmp(p1,b)<=0);
killvalue(ep);
return gerepile(av,tetpil,y);
}
GEN vecteur(ep,nmax,ch)
entree *ep;
GEN nmax;
char *ch;
{
GEN y,p1,t;
long i,m;
if((typ(nmax)!=1) || (signe(nmax)<0)) err(vecer1);
m=itos(nmax);
y=cgetg(m+1 ,17);
newvalue(ep, gun); p1 = (GEN)ep->value;
for(i=1;i<=m;i++)
{
affsi(i,p1);
t=lisexpr(ch);
y[i] = isonstack(t) ? (long)t : lcopy(t);
}
killvalue(ep);
return y;
}
GEN vvecteur(ep,nmax,ch)
entree *ep;
GEN nmax;
char *ch;
{
GEN y=vecteur(ep,nmax,ch);
settyp(y,18);
return y;
}
GEN matrice(ep1,ep2,nlig,ncol,ch)
entree *ep1,*ep2;
GEN nlig,ncol;
char *ch;
{
GEN y,z,t,p1,p2;
long i,j,m,n;
if((typ(nlig)!=1) || (signe(nlig)<=0)) err(mater1);
if((typ(ncol)!=1) || (signe(ncol)<=0)) err(mater1);
n=itos(nlig); m=itos(ncol);
newvalue(ep1, gun); p1 = (GEN)ep1->value;
newvalue(ep2, gun); p2 = (GEN)ep2->value;
y=cgetg(m+1 ,19);
for(i=1;i<=m;i++)
{
affsi(i,p2);
z=cgetg(n+1 ,18);
y[i]=(long)z;
for(j=1;j<=n;j++)
{
affsi(j,p1);
t=lisexpr(ch);
z[j] = isonstack(t) ? (long)t : lcopy(t);
}
}
killvalue(ep1); killvalue(ep2);
return y;
}
GEN divsomme(ep,num,ch)
entree *ep;
GEN num;
char *ch;
{
GEN y,z,p1;
long av=avma,tetpil,d,n,d2;
newvalue(ep, gun); p1 = (GEN)ep->value;
n=itos(num);/* provisoire */
tetpil=av=avma;y=gzero;
for(d = d2 = 1; d2 < n; d++, d2 += d+d-1)
if (!(n%d))
{
affsi(d,p1);y=gadd(y, lisexpr(ch));
affsi(n/d,p1); z = lisexpr(ch);
tetpil=avma; y=gadd(y,z);
}
if (d2 == n)
{
affsi(d,p1); z = lisexpr(ch);
tetpil=avma; y=gadd(y,z);
}
killvalue(ep);
return gerepile(av,tetpil,y);
}
/********************************************************************/
/********************************************************************/
/** **/
/** CALCUL D'UNE INTEGRALE **/
/** ( Methode de ROMBERG ) **/
/** **/
/********************************************************************/
/********************************************************************/
GEN qromb(ep,a,b,ch,prec)
entree *ep;
GEN a,b;
char *ch;
long prec;
#define JMAX 25
#define JMAXP JMAX+3
#define KLOC 5
{
GEN ss,dss,s,h,q1,p1,p2,p3,p4,p,qlint,del,sz,x,sum;
long av,tetpil,j,j1,j2,lim,l,it,tnm,sig;
av=avma;l=prec;
if(typ(a)!=2) { p=cgetr(prec);gaffect(a,p);a=p;}
if(typ(b)!=2) { p=cgetr(prec);gaffect(b,p);b=p;}
qlint=subrr(b,a);sig=signe(qlint);
if(!sig) {avma=av;return gzero;}
newvalue(ep,cgetr(prec)); q1=(GEN)ep->value;
if(sig<0) {setsigne(qlint,1);s=a;a=b;b=s;}
p3=cgetg(KLOC+1,17);p4=cgetg(KLOC+1,17);
s=cgetg(JMAXP,17);h=cgetg(JMAXP,17);
affsr(1,h[1]=lgetr(l));
gaffect(a,q1); p1=lisexpr(ch); if (!isonstack(p1)) p1=gcopy(p1);
gaffect(b,q1); p2=lisexpr(ch);
s[1]=lmul2n(gmul(qlint,gadd(p1,p2)),-1);it=1;sz=gmul(gzero,s[1]);
s[2]=s[1];h[2]=lshiftr(h[1],-2);
for(j=2;j<=JMAX;j++)
{
tnm=it;del=divrs(qlint,tnm);x=addrr(a,shiftr(del,-1));
for(sum=sz,j1=1;j1<=it;j1++,x=addrr(x,del))
{
affrr(x,q1);p1=lisexpr(ch);sum=gadd(sum,p1);
}
it*=2;
s[j]=lmul2n(gadd(s[j-1],gmul(sum,del)),-1);
if(j>=KLOC)
{
for(j1=1;j1<=KLOC;j1++)
{
p3[j1]=s[j1+j-KLOC];p4[j1]=h[j1+j-KLOC];
}
tetpil=avma;ss=polint(p4,p3,gzero,&dss);
j1=gexpo(ss);j2=gexpo(dss);lim=32*(prec-2)-j-5;
if((((j1-j2)>lim))||((j1< -lim)&&(j2<j1-1)))
{
if(gcmp0(gimag(ss))) ss=greal(ss);
tetpil=avma;
killvalue(ep);
return gerepile(av,tetpil,gmulsg(sig,ss));
}
}
s[j+1]=s[j];h[j+1]=lshiftr(h[j],-2);
}
err(intger2);
}
GEN qromo(ep,a,b,ch,prec)
entree *ep;
GEN a,b;
char *ch;
long prec;
#undef JMAX
#define JMAX 16
#define JMAXP JMAX+3
#define KLOC 5
{
GEN ss,dss,s,h,q1,sz,p1,p3,p4,p,qlint,del,ddel,x,sum;
long av,tetpil,j,j1,j2,lim,l,it,tnm,sig;
av=avma;l=prec;
if(typ(a)!=2) { p=cgetr(prec);gaffect(a,p);a=p;}
if(typ(b)!=2) { p=cgetr(prec);gaffect(b,p);b=p;}
qlint=subrr(b,a);sig=signe(qlint);
if(!sig) {avma=av;return gzero;}
if(sig<0) {setsigne(qlint,1);s=a;a=b;b=s;}
newvalue(ep,cgetr(prec)); q1=(GEN)ep->value;
p3=cgetg(KLOC+1,17);p4=cgetg(KLOC+1,17);
s=cgetg(JMAXP,17);h=cgetg(JMAXP,17);
affsr(1,h[1]=lgetr(l));affrr(shiftr(addrr(a,b),-1),q1);p1=lisexpr(ch);
s[1]=lmul(qlint,p1);it=1;sz=gmul(gzero,s[1]);
s[2]=s[1];h[2]=ldivrs(h[1],9);
for(j=2;j<=JMAX;j++)
{
tnm=it;del=divrs(qlint,3*tnm);ddel=shiftr(del,1);
x=addrr(a,shiftr(del,-1));sum=sz;
for(j1=1;j1<=it;j1++)
{
affrr(x,q1);p1=lisexpr(ch);sum=gadd(sum,p1);x=addrr(x,ddel);
affrr(x,q1);p1=lisexpr(ch);sum=gadd(sum,p1);x=addrr(x,del);
}
it*=3;
s[j]=ladd(gdivgs(s[j-1],3),gmul(sum,del));
if(j>=KLOC)
{
for(j1=1;j1<=KLOC;j1++)
{
p3[j1]=s[j1+j-KLOC];p4[j1]=h[j1+j-KLOC];
}
tetpil=avma;ss=polint(p4,p3,gzero,&dss);
j1=gexpo(ss);j2=gexpo(dss);lim=32*(prec-2)-(3*j/2)-5;
if((((j1-j2)>lim))||((j1< -lim)&&(j2<j1-1)))
{
if(gcmp0(gimag(ss))) ss=greal(ss);
tetpil=avma; killvalue(ep);
return gerepile(av,tetpil,gmulsg(sig,ss));
}
}
s[j+1]=s[j];h[j+1]=ldivrs(h[j],9);
}
err(intger2);
}
GEN qromi(ep,a,b,ch,prec)
entree *ep;
GEN a,b;
char *ch;
long prec;
#undef JMAX
#define JMAX 16
#define JMAXP JMAX+3
#define KLOC 5
{
GEN ss,dss,s,h,q1,sz,p1,p3,p4,p,qlint,del,ddel,x,sum;
long av,tetpil,j,j1,j2,lim,l,it,tnm,sig;
av=avma;l=prec;
p=cgetr(prec);gaffect(ginv(a),p);a=p;
p=cgetr(prec);gaffect(ginv(b),p);b=p;
qlint=subrr(b,a);sig= -signe(qlint);
if(!sig) {avma=av;return gzero;}
if(sig>0) {setsigne(qlint,1);s=a;a=b;b=s;}
newvalue(ep,cgetr(prec)); q1=(GEN)ep->value;
p3=cgetg(KLOC+1,17);p4=cgetg(KLOC+1,17);
s=cgetg(JMAXP,17);h=cgetg(JMAXP,17);
affsr(1,h[1]=lgetr(l));x=divsr(2,addrr(a,b));
affrr(x,q1);
p1=gmul(lisexpr(ch),mulrr(x,x));s[1]=lmul(qlint,p1);it=1;
sz=gmul(gzero,s[1]);s[2]=s[1];h[2]=ldivrs(h[1],9);
for(j=2;j<=JMAX;j++)
{
tnm=it;del=divrs(qlint,3*tnm);ddel=shiftr(del,1);
x=addrr(a,shiftr(del,-1));sum=sz;
for(j1=1;j1<=it;j1++)
{
divsrz(1,x,q1);p1=gmul(lisexpr(ch),mulrr(q1,q1));
sum=gadd(sum,p1);x=addrr(x,ddel);
divsrz(1,x,q1);p1=gmul(lisexpr(ch),mulrr(q1,q1));
sum=gadd(sum,p1);x=addrr(x,del);
}
it*=3;
s[j]=ladd(gdivgs(s[j-1],3),gmul(sum,del));
if(j>=KLOC)
{
for(j1=1;j1<=KLOC;j1++)
{
p3[j1]=s[j1+j-KLOC];p4[j1]=h[j1+j-KLOC];
}
tetpil=avma;ss=polint(p4,p3,gzero,&dss);
j1=gexpo(ss);j2=gexpo(dss);lim=32*(prec-2)-(3*j/2)-5;
if((((j1-j2)>lim))||((j1< -lim)&&(j2<j1-1)))
{
if(gcmp0(gimag(ss))) ss=greal(ss);
tetpil=avma; killvalue(ep);
return gerepile(av,tetpil,gmulsg(sig,ss));
}
}
s[j+1]=s[j];h[j+1]=ldivrs(h[j],9);
}
err(intger2);
}
GEN rombint(ep,a,b,ch,prec)
entree *ep;
GEN a,b;
char *ch;
long prec;
{
GEN s,p1,p2,p3;
static long p4[3]={0x1ff0003,0xff000003,1};
long l,av,tetpil;
l=gcmp(b,a);
if(!l) return gzero;
if(l<0) {s=a;a=b;b=s;}
av=avma;
if(gcmpgs(b,100)>=0)
{
if(gcmpgs(a,1)>=0) return qromi(ep,a,b,ch,prec);
p1=qromi(ep,gun,b,ch,prec);
if(gcmpgs(a,-100)>=0)
{
p2=qromo(ep,a,gun,ch,prec);tetpil=avma;
return gerepile(av,tetpil,gadd(p1,p2));
}
p2=qromo(ep,p4,gun,ch,prec);p3=gadd(p2,qromi(ep,a,p4,ch,prec));
tetpil=avma;return gerepile(av,tetpil,gadd(p1,p3));
}
if(gcmpgs(a,-100)>=0) return qromo(ep,a,b,ch,prec);
if(gcmpgs(b,-1)>=0)
{
p1=qromi(ep,a,p4,ch,prec);p2=qromo(ep,p4,b,ch,prec);tetpil=avma;
return gerepile(av,tetpil,gadd(p1,p2));
}
return qromi(ep,a,b,ch,prec);
}
/********************************************************************/
/********************************************************************/
/** **/
/** SOMMATION DE SERIES **/
/** **/
/********************************************************************/
/********************************************************************/
void eulsum(sum,term,jterm,tab,dsum,prec)
GEN *sum,term,tab[];
long jterm,prec,*dsum;
{
long j;
static long nterm;
GEN tmp,dum,p1;
static GEN unreel;
if(jterm==1)
{
nterm=1;tab[1]=term;*sum=gmul2n(term,-1);affsr(1,unreel=cgetr(prec));
}
else
{
tmp=tab[1];tab[1]=gmul(unreel,term);
for(j=1;j<nterm;j++)
{
dum=tab[j+1];tab[j+1]=gmul2n(gadd(tab[j],tmp),-1);tmp=dum;
}
tab[nterm+1]=gmul2n(gadd(tab[nterm],tmp),-1);
if(gcmp(gabs(tab[nterm+1],prec),gabs(tab[nterm],prec))<=0)
p1=gmul2n(tab[++nterm],-1);
else p1=tab[nterm+1];
*sum=gadd(*sum,p1);*dsum=gexpo(p1);
}
}
GEN sumalt(ep,a,ch,prec)
entree *ep;
GEN a;
char *ch;
#define JMIT 10000
{
long av,tetpil,j,nterm,jterm;
GEN p1,sum,sum0,q1,tmp,dum,*tab,unreel;
newvalue(ep,gun); q1=(GEN)ep->value; gaffect(a,q1);
tab=(GEN *)newbloc(JMIT);
av=avma;
p1=lisexpr(ch);affsr(1,unreel=cgetr(prec));
nterm=1;tab[1]=p1;sum=gmul2n(p1,-1);
for(jterm=2;jterm<=JMIT;jterm++)
{
tmp=tab[1];a=gaddsg(1,a);gaffect(a,q1);tab[1]=gmul(unreel,lisexpr(ch));
for(j=1;j<nterm;j++)
{
dum=tab[j+1];tab[j+1]=gmul2n(gadd(tab[j],tmp),-1);tmp=dum;
}
tab[nterm+1]=gmul2n(gadd(tab[nterm],tmp),-1);
if(gcmp(gabs(tab[nterm+1],prec),gabs(tab[nterm],prec))<=0)
p1=gmul2n(tab[++nterm],-1);
else p1=tab[nterm+1];
sum0=sum;sum=gadd(sum,p1);
if((gcmp0(p1)||(gexpo(p1)<-32*(prec-2)+5)||(gegal(sum,sum0)))&&(jterm>=10))
{
tetpil=avma;
killbloc((GEN)tab);
killvalue(ep);
return gerepile(av,tetpil,gcopy(sum));
}
}
err(eulsumer1);
}
GEN sumpos(ep,a,ch,prec)
entree *ep;
GEN a;
char *ch;
{
long av,tetpil,k,jterm,dsum;
GEN sum,term,q1,p1,*tab,unreel,r;
tab=(GEN *)newbloc(JMIT);
av = avma; a=gsubgs(a,1);
affsr(1,unreel=cgetr(prec));term=gzero;r=gun;k=0;
p1=cgeti(prec+1);p1[1]=0x1000001+prec;
newvalue(ep,p1); q1=(GEN)ep->value;
do
{
gaddz(r,a,q1);p1=gmul2n(gmul(unreel,lisexpr(ch)),k);
term=gadd(term,p1);r=shifti(r,1);k++;
}
while(gexpo(p1)>=(-32*(prec-2)+5));
sum=gzero;
eulsum(&sum,term,1,tab,&dsum,prec);
for(jterm=2;jterm<=JMIT;jterm++)
{
term=gzero;r=stoi(jterm);k=0;
do
{
gaddz(r,a,q1);p1=gmul2n(gmul(unreel,lisexpr(ch)),k);
term=gadd(term,p1);r=shifti(r,1);k++;
}
while(gexpo(p1)>=(-32*(prec-2)+5));
if(!odd(jterm)) term=gneg(term);
eulsum(&sum,term,jterm,tab,&dsum,prec);
if(dsum< -32*(prec-2)+5)
{
tetpil=avma;
killbloc((GEN)tab);
killvalue(ep);
return gerepile(av,tetpil,gcopy(sum));
}
}
killbloc((GEN)tab);
err(eulsumer1);
}
/********************************************************************/
/********************************************************************/
/** **/
/** TRACE GROSSIER **/
/** **/
/********************************************************************/
/********************************************************************/
GEN plot(ep,a,b,ch)
entree *ep;
GEN a,b;
char *ch;
#define ISCR 60
#define JSCR 21
#define BLANK ' '
#define ZERO '-'
#define YY '|'
#define XX '-'
#define FF 'x'
{
long av,av2,jz,j,i,sig;
GEN p1,p2,ysml,ybig,x,diff,dyj,dx,y[ISCR+1];
char scr[ISCR+1][JSCR+1], thestring[100];
sig=gcmp(b,a);if(!sig) return gnil;
av=avma; pariputc('\n');
if(sig<0) {x=a;a=b;b=x;}
newvalue(ep,cgetr(3)); x=(GEN)ep->value;
for(i=1;i<=ISCR;i++) y[i]=cgetr(3);
gaffect(a,x);
dx=gdivgs(gsub(b,a),(ISCR-1));ysml=gzero;ybig=gzero;
for(j=1;j<=JSCR;j++) scr[1][j]=scr[ISCR][j]=YY;
for(i=2;i<ISCR;i++)
{
scr[i][1]=scr[i][JSCR]=XX;
for(j=2;j<JSCR;j++) scr[i][j]=BLANK;
}
av2=avma;
for(i=1;i<=ISCR;i++)
{
gaffect(lisexpr(ch),y[i]);
if(gcmp(y[i],ysml)<0) ysml=y[i];
if(gcmp(y[i],ybig)>0) ybig=y[i];
gaddz(x,dx,x);avma=av2;
}
diff=gsub(ybig,ysml);
if(gcmp0(diff)) {ybig=gaddsg(1,ybig);diff=gun;}
dyj=gdivsg(JSCR-1,diff);jz=1-itos(ground(gmul(ysml,dyj)));
av2=avma;
for(i=1;i<=ISCR;i++)
{
scr[i][jz]=ZERO;j=1+itos(ground(gmul(gsub(y[i],ysml),dyj)));
scr[i][j]=FF;avma=av2;
}
p1=cgetr(4);gaffect(ybig,p1);
sprintf(thestring, " %10.3lf ",rtodbl(p1)); pariputs(thestring);
for(i=1;i<=ISCR;i++) pariputc(scr[i][JSCR]); pariputc('\n');
for(j=(JSCR-1);j>=2;j--)
{
pariputs(" ");
for(i=1;i<=ISCR;i++) pariputc(scr[i][j]);
pariputc('\n');
}
p1=cgetr(4);gaffect(ysml,p1);
sprintf(thestring, " %10.3lf ",rtodbl(p1)); pariputs(thestring);
for(i=1;i<=ISCR;i++) pariputc(scr[i][1]);
pariputc('\n');
p1=cgetr(4);p2=cgetr(4);gaffect(a,p1);gaffect(b,p2);
sprintf(thestring, "%8s %10.3lf %44s %10.3lf\n"," ",rtodbl(p1)," ",rtodbl(p2));
pariputs(thestring);
killvalue(ep);
avma=av; pariputc('\n'); return gnil;
}
/********************************************************************/
/********************************************************************/
/** **/
/** RECHERCHE DE ZEROS REELS **/
/** **/
/********************************************************************/
/********************************************************************/
GEN zbrent(ep,a,b,ch,prec)
entree *ep;
GEN a,b;
char *ch;
long prec;
{
long av=avma,tetpil,sig,iter,itmax;
GEN q1,c,d,e,tol,toli,min1,min2,fa,fb,fc,p,q,r,s,xm;
if(typ(a)!=2) {p=cgetr(prec);gaffect(a,p);a=p;}
if(typ(b)!=2) {p=cgetr(prec);gaffect(b,p);b=p;}
sig=cmprr(b,a);if(!sig) {avma = av; return gzero;}
if(sig<0) {c=a;a=b;b=c;}
newvalue(ep,a); fa=lisexpr(ch);
changevalue(ep,b); q1=(GEN)ep->value; fb=lisexpr(ch);
if(gsigne(fa)*gsigne(fb)>0) err(zbrenter1);
itmax=32*prec+50;affsr(1,tol=cgetr(3));tol=shiftr(tol,-(32*(prec-2)-5));
fc=fb;
for(iter=1;iter<=itmax;iter++)
{
if(gsigne(fb)*gsigne(fc)>0)
{
c=a;fc=fa;d=subrr(b,a);e=d;
}
if(gcmp(gabs(fc),gabs(fb))<0)
{
a=b;b=c;c=a;fa=fb;fb=fc;fc=fa;
}
toli=mulrr(tol,absr(b));
xm=shiftr(subrr(c,b),-1);
if((cmprr(absr(xm),toli)<=0)||gcmp0(fb))
{
tetpil=avma;
killvalue(ep);
return gerepile(av,tetpil,gcopy(b));
}
if((cmprr(absr(e),toli)>=0)&&(gcmp(gabs(fa),gabs(fb))>0))
{
s=gdiv(fb,fa);
if(cmprr(a,b)==0)
{
p=gmul2n(gmul(xm,s),1);q=gsubsg(1,s);
}
else
{
q=gdiv(fa,fc);r=gdiv(fb,fc);
p=gmul2n(gmul(gsub(q,r),gmul(xm,q)),1);
p=gmul(s,gsub(p,gmul(gsub(b,a),gsubgs(r,1))));
q=gmul(gmul(gsubgs(q,1),gsubgs(r,1)),gsubgs(s,1));
}
if(gsigne(p)>0) q=gneg(q);
p=gabs(p);
min1=gsub(gmulsg(3,gmul(xm,q)),gabs(gmul(q,toli)));
min2=gabs(gmul(e,q));
if(gcmp(gmul2n(p,1),gmin(min1,min2))<0) { e=d;d=gdiv(p,q);}
else {d=xm;e=d;}
}
else {d=xm;e=d;}
a=b;fa=fb;
if(gcmp(gabs(d),toli)>0) b=addrr(b,d);
else
{
if(gsigne(xm)>0) b=addrr(b,absr(toli));
else b=subrr(b,absr(toli));
}
gaffect(b,q1); ;fb=lisexpr(ch);
}
err(zbrenter2);
}